home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C / Utilities / Calc 1.24.7 / lib / pell.cal < prev    next >
Encoding:
Text File  |  1992-02-24  |  1.3 KB  |  74 lines  |  [TEXT/R*ch]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
  7.  * Type the solution to pells equation for a particular D.
  8.  */
  9.  
  10. define pell(D)
  11. {
  12.     local X, Y;
  13.  
  14.     X = pellx(D);
  15.     if (isnull(X)) {
  16.         print "D=":D:" is square";
  17.         return;
  18.     }
  19.     Y = isqrt((X^2 - 1) / D);
  20.     print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2;
  21. }
  22.  
  23.  
  24. /*
  25.  * Function to solve Pell's equation
  26.  * Returns the solution X to:
  27.  *    X^2 - D * Y^2 = 1
  28.  */
  29. define pellx(D)
  30. {
  31.     local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n, ans, tmp;
  32.  
  33.     mat ans[2,2];
  34.     mat tmp[2,2];
  35.  
  36.     R = isqrt(D);
  37.     Vp = D - R^2;
  38.     if (Vp == 0)
  39.         return;
  40.     Rp = R + R;
  41.     U = Rp;
  42.     Up = U;
  43.     V = 1;
  44.     A = 0;
  45.     n = 0;
  46.     ans[0,0] = 1;
  47.     ans[1,1] = 1;
  48.     tmp[0,1] = 1;
  49.     tmp[1,0] = 1;
  50.     do {
  51.         T = V;
  52.         V = A * (Up - U) + Vp;
  53.         Vp = T;
  54.         A = U // V;
  55.         Up = U;
  56.         U = Rp - U % V;
  57.         tmp[0,0] = A;
  58.         ans *= tmp;
  59.         n++;
  60.     } while (A != Rp);
  61.     Q2 = ans[[1]];
  62.     Q1 = isqrt(Q2^2 * D + 1);
  63.     if (isodd(n)) {
  64.         T = Q1^2 + D * Q2^2;
  65.         Q2 = Q1 * Q2 * 2;
  66.         Q1 = T;
  67.     }
  68.     return Q1;
  69. }
  70.  
  71. global lib_debug;
  72. if (!isnum(lib_debug) || lib_debug>0) print "pell(D) defined";
  73. if (!isnum(lib_debug) || lib_debug>0) print "pellx(D) defined";
  74.